function vmr = mergehemivmrs(h1, h2)
% mergehemivmrs  - merge segmented VMRs of LH and RH
%
% FORMAT:       vmr = mergehemivmrs(h1, h2)
%
% Input fields:
%
%       h1, h2      segmented hemisphere VMR
%
% Output fields:
%
%       vmr         combined VMR

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% copy object
vmr = h1.CopyObject;

% get voxels to copy
copyb = conv3d(vmr.VMRData == 0, 1);

% copy data
vmr.VMRData(copyb) = h2.VMRData(copyb);

% find hemisphere cut
po = squeeze(sum(sum(vmr.VMRData(:, :, 113:144) >= 235, 1), 2));
[pm{1:2}] = min(po);
pm = pm{2} + 112;

% get adjacent slices
h1s = vmr.VMRData(:, :, pm - 1);
h2s = vmr.VMRData(:, :, pm + 1);

% get biggest "chunk"
[cc{1:2}] = clustercoords(cat(3, h1s, h2s) > 0, 'edge');
if isempty(cc{2})
    error( ...
        'BVQXtools:InternalError', ...
        'Error generating cluster list in connected hemisphere.' ...
    );
end
cc = cc{2};
cn = 0;
ci = 0;
for c = 1:numel(cc)
    if cn < size(cc{c}, 1)
        ci = c;
        cn = size(cc{c}, 1);
    end
end

% reconnect by marking overlap
cc = cc{ci};
s1i = cc(:, 3) == 1;
mi = intersect( ...
    sub2ind(size(h1s), cc(s1i, 1), cc(s1i, 2)), ...
    sub2ind(size(h1s), cc(~s1i, 1), cc(~s1i, 2))) + ...
    (pm - 1) * numel(h1s);
vmr.VMRData(mi) = 240;

% prepare again
bless(vmr.PrepareForReco, 1);
